!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine mlw (ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy,                                                                    &
   nx, ny, ng, n_neighs, dx, dy, dt, gm, cd, cfl, div_clean_flag,                                                                  &
   fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd,                                       &
   fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny, ng, n_neighs, div_clean_flag
   real(8), intent(in) :: dx, dy, dt, gm, cd, cfl
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, az, po, gx, gy
   real(8), dimension(ny, n_neighs), intent(inout) ::                                                                              &
      fxro_bnd, fxmx_bnd, fxmy_bnd, fxmz_bnd, fxbx_bnd, fxby_bnd, fxbz_bnd, fxen_bnd, fxpo_bnd
   real(8), dimension(nx, n_neighs), intent(inout) ::                                                                              &
      fyro_bnd, fymx_bnd, fymy_bnd, fymz_bnd, fybx_bnd, fyby_bnd, fybz_bnd, fyen_bnd, fypo_bnd

   integer(4) :: i, j, ip1, jp1, im1, jm1, nxmng, nymng

   real(8) :: eps, dt_dx, dt_dy, gmm1, dx2, dy2, ch, ch2, cp2, pr, vx, vy, vz, b2, v2, art_amp, dvx1, dvx2, dvy1, dvy2
   real(8) :: div_b, s_mx, s_my, s_mz, s_bx, s_by, s_bz, s_en

   real(8), dimension(nx, ny) :: ro_pre, mx_pre, my_pre, mz_pre, bx_pre, by_pre, bz_pre, en_pre, po_pre, az_pre
   real(8), dimension(nx, ny) :: fxro1, fxmx1, fxmy1, fxmz1, fxbx1, fxby1, fxbz1, fxen1, fxpo1
   real(8), dimension(nx, ny) :: fyro1, fymx1, fymy1, fymz1, fybx1, fyby1, fybz1, fyen1, fypo1
   real(8), dimension(nx, ny) :: fxro2, fxmx2, fxmy2, fxmz2, fxbx2, fxby2, fxbz2, fxen2, fxpo2
   real(8), dimension(nx, ny) :: fyro2, fymx2, fymy2, fymz2, fybx2, fyby2, fybz2, fyen2, fypo2

!----------------------------------------------------------------------------------------------------------------------------------|
   eps = 1.0d-12
   dt_dx = dt / dx
   dt_dy = dt / dy
   gmm1 = gm - 1.0d0
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   ch = cfl * min (dx, dy) / dt * div_clean_flag
   ch2 = ch * ch
!   cp2 = ch * cd / 6.0d0 + eps
   cp2 = -dt * ch2 / log (cd) + eps

   do j = 1, ny
   do i = 1, nx
      vx = mx(i, j) / ro(i, j)
      vy = my(i, j) / ro(i, j)
      vz = mz(i, j) / ro(i, j)
      b2 = bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en(i, j) - v2 * ro(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1

      fxro1(i, j) = mx(i, j)
	   fxmx1(i, j) = vx * mx(i, j) + pr + b2 / 2.0d0 - bx(i, j) * bx(i, j)
	   fxmy1(i, j) = vx * my(i, j) - bx(i, j) * by(i, j)
	   fxmz1(i, j) = vx * mz(i, j) - bx(i, j) * bz(i, j)
	   fxbx1(i, j) = po(i, j)
	   fxby1(i, j) = vx * by(i, j) - bx(i, j) * vy
	   fxbz1(i, j) = vx * bz(i, j) - bx(i, j) * vz
	   fxen1(i, j) = (en(i, j) + pr + b2 / 2.0d0) * vx - (bx(i, j) * vx + by(i, j) * vy + bz(i, j) * vz) * bx(i, j)
      fxpo1(i, j) = ch2 * bx(i, j)

      fyro1(i, j) = my(i, j)
	   fymx1(i, j) = vy * mx(i, j) - by(i, j) * bx(i, j)
	   fymy1(i, j) = vy * my(i, j) + pr + b2 / 2.0d0 - by(i, j) * by(i, j)
	   fymz1(i, j) = vy * mz(i, j) - by(i, j) * bz(i, j)
	   fybx1(i, j) = vy * bx(i, j) - by(i, j) * vx
	   fyby1(i, j) = po(i, j)
	   fybz1(i, j) = vy * bz(i, j) - by(i, j) * vz
	   fyen1(i, j) = (en(i, j) + pr + b2 / 2.0d0) * vy - (bx(i, j) * vx + by(i, j) * vy + bz(i, j) * vz) * by(i, j)
      fypo1(i, j) = ch2 * by(i, j)
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  predictor-corrector step
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny - 1
   do i = 1, nx - 1
      ip1 = i + 1
      jp1 = j + 1
      ro_pre(i, j) = (ro(i, j) + ro(ip1, j) + ro(ip1, jp1) + ro(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxro1(ip1, j) + fxro1(ip1, jp1) - fxro1(i, j) - fxro1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fyro1(i, jp1) + fyro1(ip1, jp1) - fyro1(i, j) - fyro1(ip1, j)) / 2.0d0
      mx_pre(i, j) = (mx(i, j) + mx(ip1, j) + mx(ip1, jp1) + mx(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxmx1(ip1, j) + fxmx1(ip1, jp1) - fxmx1(i, j) - fxmx1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fymx1(i, jp1) + fymx1(ip1, jp1) - fymx1(i, j) - fymx1(ip1, j)) / 2.0d0
      my_pre(i, j) = (my(i, j) + my(ip1, j) + my(ip1, jp1) + my(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxmy1(ip1, j) + fxmy1(ip1, jp1) - fxmy1(i, j) - fxmy1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fymy1(i, jp1) + fymy1(ip1, jp1) - fymy1(i, j) - fymy1(ip1, j)) / 2.0d0
      mz_pre(i, j) = (mz(i, j) + mz(ip1, j) + mz(ip1, jp1) + mz(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxmz1(ip1, j) + fxmz1(ip1, jp1) - fxmz1(i, j) - fxmz1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fymz1(i, jp1) + fymz1(ip1, jp1) - fymz1(i, j) - fymz1(ip1, j)) / 2.0d0
      bx_pre(i, j) = (bx(i, j) + bx(ip1, j) + bx(ip1, jp1) + bx(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxbx1(ip1, j) + fxbx1(ip1, jp1) - fxbx1(i, j) - fxbx1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fybx1(i, jp1) + fybx1(ip1, jp1) - fybx1(i, j) - fybx1(ip1, j)) / 2.0d0
      by_pre(i, j) = (by(i, j) + by(ip1, j) + by(ip1, jp1) + by(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxby1(ip1, j) + fxby1(ip1, jp1) - fxby1(i, j) - fxby1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fyby1(i, jp1) + fyby1(ip1, jp1) - fyby1(i, j) - fyby1(ip1, j)) / 2.0d0
      bz_pre(i, j) = (bz(i, j) + bz(ip1, j) + bz(ip1, jp1) + bz(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxbz1(ip1, j) + fxbz1(ip1, jp1) - fxbz1(i, j) - fxbz1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fybz1(i, jp1) + fybz1(ip1, jp1) - fybz1(i, j) - fybz1(ip1, j)) / 2.0d0
      en_pre(i, j) = (en(i, j) + en(ip1, j) + en(ip1, jp1) + en(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxen1(ip1, j) + fxen1(ip1, jp1) - fxen1(i, j) - fxen1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fyen1(i, jp1) + fyen1(ip1, jp1) - fyen1(i, j) - fyen1(ip1, j)) / 2.0d0
      po_pre(i, j) = (po(i, j) + po(ip1, j) + po(ip1, jp1) + po(i, jp1)) / 4.0d0 -                                                 &
         dt_dx * (fxpo1(ip1, j) + fxpo1(ip1, jp1) - fxpo1(i, j) - fxpo1(i, jp1)) / 2.0d0 -                                         &
         dt_dy * (fypo1(i, jp1) + fypo1(ip1, jp1) - fypo1(i, j) - fypo1(ip1, j)) / 2.0d0
   enddo
   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      div_b = (bx(i + 1, j) - bx(i - 1, j)) / dx2 * div_clean_flag + (by(i, j + 1) - by(i, j - 1)) / dy2 * div_clean_flag
      s_mx = -div_b * bx(i, j) + ro(i, j) * gx(i, j)
      s_my = -div_b * by(i, j) + ro(i, j) * gy(i, j)
      s_mz = -div_b * bz(i, j)
      s_bx = -div_b * mx(i, j) / ro(i, j)
      s_by = -div_b * my(i, j) / ro(i, j)
      s_bz = -div_b * mz(i, j) / ro(i, j)
      s_en = -div_b * (bx(i, j) * mx(i, j) + by(i, j) * my(i, j) + bz(i, j) * mz(i, j)) / ro(i, j) -                               &
              (po(i + 1, j) - po(i - 1, j)) / dx2 * bx(i, j) * div_clean_flag -                                                    &
              (po(i, j + 1) - po(i, j - 1)) / dy2 * by(i, j) * div_clean_flag +                                                    &
              mx(i, j) * gx(i, j) + my(i, j) * gy(i, j)
!      s_en = -(po(i + 1, j) - po(i - 1, j)) / dx2 * bx(i, j) * div_clean_flag -                                                    &
!              (po(i, j + 1) - po(i, j - 1)) / dy2 * by(i, j) * div_clean_flag +                                                    &
!              mx(i, j) * gx(i, j) + my(i, j) * gy(i, j)

      mx_pre(i, j) = mx_pre(i, j) + s_mx * dt * 0.5d0
      my_pre(i, j) = my_pre(i, j) + s_my * dt * 0.5d0
      mz_pre(i, j) = mz_pre(i, j) + s_mz * dt * 0.5d0
      bx_pre(i, j) = bx_pre(i, j) + s_bx * dt * 0.5d0
      by_pre(i, j) = by_pre(i, j) + s_by * dt * 0.5d0
      bz_pre(i, j) = bz_pre(i, j) + s_bz * dt * 0.5d0
      en_pre(i, j) = en_pre(i, j) + s_en * dt * 0.5d0
      po_pre(i, j) = po_pre(i, j) * exp (-dt * 0.5d0 * ch2 / cp2 * div_clean_flag)
      az_pre(i, j) = az(i, j) + (mx(i, j) * by(i, j) - my(i, j) * bx(i, j)) / ro(i, j) * dt * 0.5d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  final result
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      vx = mx_pre(i, j) / ro_pre(i, j)
      vy = my_pre(i, j) / ro_pre(i, j)
      vz = mz_pre(i, j) / ro_pre(i, j)
      b2 = bx_pre(i, j) * bx_pre(i, j) + by_pre(i, j) * by_pre(i, j) + bz_pre(i, j) * bz_pre(i, j)
      v2 = vx * vx + vy * vy + vz * vz
      pr = (en_pre(i, j) - v2 * ro_pre(i, j) / 2.0d0 - b2 / 2.0d0) * gmm1

      fxro2(i, j) = mx_pre(i, j)
	   fxmx2(i, j) = vx * mx_pre(i, j) + pr + b2 / 2.0d0 - bx_pre(i, j) * bx_pre(i, j)
	   fxmy2(i, j) = vx * my_pre(i, j) - bx_pre(i, j) * by_pre(i, j)
	   fxmz2(i, j) = vx * mz_pre(i, j) - bx_pre(i, j) * bz_pre(i, j)
	   fxbx2(i, j) = po_pre(i, j)
	   fxby2(i, j) = vx * by_pre(i, j) - bx_pre(i, j) * vy
	   fxbz2(i, j) = vx * bz_pre(i, j) - bx_pre(i, j) * vz
	   fxen2(i, j) = (en_pre(i, j) + pr + b2 / 2.0d0) * vx -                                                                        &
         (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * bx_pre(i, j)
      fxpo2(i, j) = ch2 * bx_pre(i, j)

      fyro2(i, j) = my_pre(i, j)
	   fymx2(i, j) = vy * mx_pre(i, j) - by_pre(i, j) * bx_pre(i, j)
	   fymy2(i, j) = vy * my_pre(i, j) + pr + b2 / 2.0d0 - by_pre(i, j) * by_pre(i, j)
	   fymz2(i, j) = vy * mz_pre(i, j) - by_pre(i, j) * bz_pre(i, j)
	   fybx2(i, j) = vy * bx_pre(i, j) - by_pre(i, j) * vx
	   fyby2(i, j) = po_pre(i, j)
	   fybz2(i, j) = vy * bz_pre(i, j) - by_pre(i, j) * vz
	   fyen2(i, j) = (en_pre(i, j) + pr + b2 / 2.0d0) * vy -                                                                        &
         (bx_pre(i, j) * vx + by_pre(i, j) * vy + bz_pre(i, j) * vz) * by_pre(i, j)
      fypo2(i, j) = ch2 * by_pre(i, j)
   enddo
   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      ip1 = i + 1
      jp1 = j + 1
      im1 = i - 1
      jm1 = j - 1
      ro_pre(i, j) = ro(i, j) -                                                                                                    &
         dt_dx * ((fxro1(ip1, j) - fxro1(im1, j)) / 2.0d0 +                                                                        &
         (fxro2(i, j) + fxro2(i, jm1) - fxro2(im1, j) - fxro2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fyro1(i, jp1) - fyro1(i, jm1)) / 2.0d0 +                                                                        &
         (fyro2(i, j) + fyro2(im1, j) - fyro2(i, jm1) - fyro2(im1, jm1)) / 2.0d0) / 2.0d0
      mx_pre(i, j) = mx(i, j) -                                                                                                    &
         dt_dx * ((fxmx1(ip1, j) - fxmx1(im1, j)) / 2.0d0 +                                                                        &
         (fxmx2(i, j) + fxmx2(i, jm1) - fxmx2(im1, j) - fxmx2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fymx1(i, jp1) - fymx1(i, jm1)) / 2.0d0 +                                                                        &
         (fymx2(i, j) + fymx2(im1, j) - fymx2(i, jm1) - fymx2(im1, jm1)) / 2.0d0) / 2.0d0
      my_pre(i, j) = my(i, j) -                                                                                                    &
         dt_dx * ((fxmy1(ip1, j) - fxmy1(im1, j)) / 2.0d0 +                                                                        &
         (fxmy2(i, j) + fxmy2(i, jm1) - fxmy2(im1, j) - fxmy2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fymy1(i, jp1) - fymy1(i, jm1)) / 2.0d0 +                                                                        &
         (fymy2(i, j) + fymy2(im1, j) - fymy2(i, jm1) - fymy2(im1, jm1)) / 2.0d0) / 2.0d0
      mz_pre(i, j) = mz(i, j) -                                                                                                    &
         dt_dx * ((fxmz1(ip1, j) - fxmz1(im1, j)) / 2.0d0 +                                                                        &
         (fxmz2(i, j) + fxmz2(i, jm1) - fxmz2(im1, j) - fxmz2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fymz1(i, jp1) - fymz1(i, jm1)) / 2.0d0 +                                                                        &
         (fymz2(i, j) + fymz2(im1, j) - fymz2(i, jm1) - fymz2(im1, jm1)) / 2.0d0) / 2.0d0
      bx_pre(i, j) = bx(i, j) -                                                                                                    &
         dt_dx * ((fxbx1(ip1, j) - fxbx1(im1, j)) / 2.0d0 +                                                                        &
         (fxbx2(i, j) + fxbx2(i, jm1) - fxbx2(im1, j) - fxbx2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fybx1(i, jp1) - fybx1(i, jm1)) / 2.0d0 +                                                                        &
         (fybx2(i, j) + fybx2(im1, j) - fybx2(i, jm1) - fybx2(im1, jm1)) / 2.0d0) / 2.0d0
      by_pre(i, j) = by(i, j) -                                                                                                    &
         dt_dx * ((fxby1(ip1, j) - fxby1(im1, j)) / 2.0d0 +                                                                        &
         (fxby2(i, j) + fxby2(i, jm1) - fxby2(im1, j) - fxby2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fyby1(i, jp1) - fyby1(i, jm1)) / 2.0d0 +                                                                        &
         (fyby2(i, j) + fyby2(im1, j) - fyby2(i, jm1) - fyby2(im1, jm1)) / 2.0d0) / 2.0d0
      bz_pre(i, j) = bz(i, j) -                                                                                                    &
         dt_dx * ((fxbz1(ip1, j) - fxbz1(im1, j)) / 2.0d0 +                                                                        &
         (fxbz2(i, j) + fxbz2(i, jm1) - fxbz2(im1, j) - fxbz2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fybz1(i, jp1) - fybz1(i, jm1)) / 2.0d0 +                                                                        &
         (fybz2(i, j) + fybz2(im1, j) - fybz2(i, jm1) - fybz2(im1, jm1)) / 2.0d0) / 2.0d0
      en_pre(i, j) = en(i, j) -                                                                                                    &
         dt_dx * ((fxen1(ip1, j) - fxen1(im1, j)) / 2.0d0 +                                                                        &
         (fxen2(i, j) + fxen2(i, jm1) - fxen2(im1, j) - fxen2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fyen1(i, jp1) - fyen1(i, jm1)) / 2.0d0 +                                                                        &
         (fyen2(i, j) + fyen2(im1, j) - fyen2(i, jm1) - fyen2(im1, jm1)) / 2.0d0) / 2.0d0
      po_pre(i, j) = po(i, j) -                                                                                                    &
         dt_dx * ((fxpo1(ip1, j) - fxpo1(im1, j)) / 2.0d0 +                                                                        &
         (fxpo2(i, j) + fxpo2(i, jm1) - fxpo2(im1, j) - fxpo2(im1, jm1)) / 2.0d0) / 2.0d0 -                                        &
         dt_dy * ((fypo1(i, jp1) - fypo1(i, jm1)) / 2.0d0 +                                                                        &
         (fypo2(i, j) + fypo2(im1, j) - fypo2(i, jm1) - fypo2(im1, jm1)) / 2.0d0) / 2.0d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  artificial viscosity
!----------------------------------------------------------------------------------------------------------------------------------|
   art_amp = 3.0d0
   do j = 2, ny - 1
   do i = 2, nx - 1
      ip1 = i + 1
      im1 = i - 1
      jp1 = j + 1
      jm1 = j - 1
      dvx1 = art_amp * abs (mx(ip1, j) / ro(ip1, j) - mx(i, j) / ro(i, j))
      dvx2 = art_amp * abs (mx(i, j) / ro(i, j) - mx(im1, j) / ro(im1, j))
      dvy1 = art_amp * abs (my(i, jp1) / ro(i, jp1) - my(i, j) / ro(i, j))
      dvy2 = art_amp * abs (my(i, j) / ro(i, j) - my(i, jm1) / ro(i, jm1))
      ro(i, j) = ro_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (ro(ip1, j) - ro(i, j)) - dvx2 * (ro(i, j) - ro(im1, j))) +                                               &
         dt_dy * (dvy1 * (ro(i, jp1) - ro(i, j)) - dvy2 * (ro(i, j) - ro(i, jm1)))
      mx(i, j) = mx_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (mx(ip1, j) - mx(i, j)) - dvx2 * (mx(i, j) - mx(im1, j))) +                                               &
         dt_dy * (dvy1 * (mx(i, jp1) - mx(i, j)) - dvy2 * (mx(i, j) - mx(i, jm1)))
      my(i, j) = my_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (my(ip1, j) - my(i, j)) - dvx2 * (my(i, j) - my(im1, j))) +                                               &
         dt_dy * (dvy1 * (my(i, jp1) - my(i, j)) - dvy2 * (my(i, j) - my(i, jm1)))
      mz(i, j) = mz_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (mz(ip1, j) - mz(i, j)) - dvx2 * (mz(i, j) - mz(im1, j))) +                                               &
         dt_dy * (dvy1 * (mz(i, jp1) - mz(i, j)) - dvy2 * (mz(i, j) - mz(i, jm1)))
      bx(i, j) = bx_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (bx(ip1, j) - bx(i, j)) - dvx2 * (bx(i, j) - bx(im1, j))) +                                               &
         dt_dy * (dvy1 * (bx(i, jp1) - bx(i, j)) - dvy2 * (bx(i, j) - bx(i, jm1)))
      by(i, j) = by_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (by(ip1, j) - by(i, j)) - dvx2 * (by(i, j) - by(im1, j))) +                                               &
         dt_dy * (dvy1 * (by(i, jp1) - by(i, j)) - dvy2 * (by(i, j) - by(i, jm1)))
      bz(i, j) = bz_pre(i, j) +                                                                                                    &
         dt_dx * (dvx1 * (bz(ip1, j) - bz(i, j)) - dvx2 * (bz(i, j) - bz(im1, j))) +                                               &
         dt_dy * (dvy1 * (bz(i, jp1) - bz(i, j)) - dvy2 * (bz(i, j) - bz(i, jm1)))
      en(i, j) = en_pre(i, j) +                                                                                                        &
         dt_dx * (dvx1 * (en(ip1, j) - en(i, j)) - dvx2 * (en(i, j) - en(im1, j))) +                                               &
         dt_dy * (dvy1 * (en(i, jp1) - en(i, j)) - dvy2 * (en(i, j) - en(i, jm1)))
      po(i, j) = po_pre(i, j) +                                                                                                        &
         dt_dx * (dvx1 * (po(ip1, j) - po(i, j)) - dvx2 * (po(i, j) - po(im1, j))) +                                               &
         dt_dy * (dvy1 * (po(i, jp1) - po(i, j)) - dvy2 * (po(i, j) - po(i, jm1)))
   enddo
   enddo

   do j = 2, ny - 1
   do i = 2, nx - 1
      div_b = (bx_pre(i + 1, j) - bx_pre(i - 1, j)) / dx2 * div_clean_flag +                                                       &
              (by_pre(i, j + 1) - by_pre(i, j - 1)) / dy2 * div_clean_flag
      s_mx = -div_b * bx_pre(i, j) + ro_pre(i, j) * gx(i, j)
      s_my = -div_b * by_pre(i, j) + ro_pre(i, j) * gy(i, j)
      s_mz = -div_b * bz_pre(i, j)
      s_bx = -div_b * mx_pre(i, j) / ro_pre(i, j)
      s_by = -div_b * my_pre(i, j) / ro_pre(i, j)
      s_bz = -div_b * mz_pre(i, j) / ro_pre(i, j)
      s_en = -div_b * (bx_pre(i, j) * mx_pre(i, j) + by_pre(i, j) * my_pre(i, j) + bz_pre(i, j) * mz_pre(i, j)) / ro_pre(i, j) -   &
              (po_pre(i + 1, j) - po_pre(i - 1, j)) / dx2 * bx_pre(i, j) * div_clean_flag -                                        &
              (po_pre(i, j + 1) - po_pre(i, j - 1)) / dy2 * by_pre(i, j) * div_clean_flag +                                        &
              mx_pre(i, j) * gx(i, j) + my_pre(i, j) * gy(i, j)
!      s_en = -(po_pre(i + 1, j) - po_pre(i - 1, j)) / dx2 * bx_pre(i, j) * div_clean_flag -                                        &
!              (po_pre(i, j + 1) - po_pre(i, j - 1)) / dy2 * by_pre(i, j) * div_clean_flag +                                        &
!              mx_pre(i, j) * gx(i, j) + my_pre(i, j) * gy(i, j)

      mx(i, j) = mx(i, j) + s_mx * dt
      my(i, j) = my(i, j) + s_my * dt
      mz(i, j) = mz(i, j) + s_mz * dt
      bx(i, j) = bx(i, j) + s_bx * dt
      by(i, j) = by(i, j) + s_by * dt
      bz(i, j) = bz(i, j) + s_bz * dt
      en(i, j) = en(i, j) + s_en * dt
      po(i, j) = po(i, j) * exp (-dt * ch2 / cp2 * div_clean_flag)
      az(i, j) = az(i, j) + (mx_pre(i, j) * by_pre(i, j) - my_pre(i, j) * bx_pre(i, j)) / ro_pre(i, j) * dt
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  bnd value
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 2, ny
      jm1 = j - 1
      nxmng = nx - ng
      fxro_bnd(j, 1) = (fxro2(ng, j) + fxro2(ng, jm1)) / 2.0d0
      fxro_bnd(j, 2) = (fxro2(nxmng, j) + fxro2(nxmng, jm1)) / 2.0d0

      fxmx_bnd(j, 1) = (fxmx2(ng, j) + fxmx2(ng, jm1)) / 2.0d0
      fxmx_bnd(j, 2) = (fxmx2(nxmng, j) + fxmx2(nxmng, jm1)) / 2.0d0

      fxmy_bnd(j, 1) = (fxmy2(ng, j) + fxmy2(ng, jm1)) / 2.0d0
      fxmy_bnd(j, 2) = (fxmy2(nxmng, j) + fxmy2(nxmng, jm1)) / 2.0d0

      fxmz_bnd(j, 1) = (fxmz2(ng, j) + fxmz2(ng, jm1)) / 2.0d0
      fxmz_bnd(j, 2) = (fxmz2(nxmng, j) + fxmz2(nxmng, jm1)) / 2.0d0

      fxbx_bnd(j, 1) = (fxbx2(ng, j) + fxbx2(ng, jm1)) / 2.0d0
      fxbx_bnd(j, 2) = (fxbx2(nxmng, j) + fxbx2(nxmng, jm1)) / 2.0d0

      fxby_bnd(j, 1) = (fxby2(ng, j) + fxby2(ng, jm1)) / 2.0d0
      fxby_bnd(j, 2) = (fxby2(nxmng, j) + fxby2(nxmng, jm1)) / 2.0d0

      fxbz_bnd(j, 1) = (fxbz2(ng, j) + fxbz2(ng, jm1)) / 2.0d0
      fxbz_bnd(j, 2) = (fxbz2(nxmng, j) + fxbz2(nxmng, jm1)) / 2.0d0

      fxen_bnd(j, 1) = (fxen2(ng, j) + fxen2(ng, jm1)) / 2.0d0
      fxen_bnd(j, 2) = (fxen2(nxmng, j) + fxen2(nxmng, jm1)) / 2.0d0

      fxpo_bnd(j, 1) = (fxpo2(ng, j) + fxpo2(ng, jm1)) / 2.0d0
      fxpo_bnd(j, 2) = (fxpo2(nxmng, j) + fxpo2(nxmng, jm1)) / 2.0d0
   enddo

   do i = 2, nx
      im1 = i - 1
      nymng = ny - ng
      fyro_bnd(i, 1) = (fyro2(i, ng) + fyro2(im1, ng)) / 2.0d0
      fyro_bnd(i, 2) = (fyro2(i, nymng) + fyro2(im1, nymng)) / 2.0d0

      fymx_bnd(i, 1) = (fymx2(i, ng) + fymx2(im1, ng)) / 2.0d0
      fymx_bnd(i, 2) = (fymx2(i, nymng) + fymx2(im1, nymng)) / 2.0d0

      fymy_bnd(i, 1) = (fymy2(i, ng) + fymy2(im1, ng)) / 2.0d0
      fymy_bnd(i, 2) = (fymy2(i, nymng) + fymy2(im1, nymng)) / 2.0d0

      fymz_bnd(i, 1) = (fymz2(i, ng) + fymz2(im1, ng)) / 2.0d0
      fymz_bnd(i, 2) = (fymz2(i, nymng) + fymz2(im1, nymng)) / 2.0d0

      fybx_bnd(i, 1) = (fybx2(i, ng) + fybx2(im1, ng)) / 2.0d0
      fybx_bnd(i, 2) = (fybx2(i, nymng) + fybx2(im1, nymng)) / 2.0d0

      fyby_bnd(i, 1) = (fyby2(i, ng) + fyby2(im1, ng)) / 2.0d0
      fyby_bnd(i, 2) = (fyby2(i, nymng) + fyby2(im1, nymng)) / 2.0d0

      fybz_bnd(i, 1) = (fybz2(i, ng) + fybz2(im1, ng)) / 2.0d0
      fybz_bnd(i, 2) = (fybz2(i, nymng) + fybz2(im1, nymng)) / 2.0d0

      fyen_bnd(i, 1) = (fyen2(i, ng) + fyen2(im1, ng)) / 2.0d0
      fyen_bnd(i, 2) = (fyen2(i, nymng) + fyen2(im1, nymng)) / 2.0d0

      fypo_bnd(i, 1) = (fypo2(i, ng) + fypo2(im1, ng)) / 2.0d0
      fypo_bnd(i, 2) = (fypo2(i, nymng) + fypo2(im1, nymng)) / 2.0d0
   enddo

   return
end subroutine mlw
